%% Navier-Stokes equation - MAC solver for structured non-regular spacing
% Documentation
%
%

%% Problem definition

addpath('source')
clear all
clc
close all
global NI
global NJ
Lx = 12.0
Ly = 2.0
h = 0.25
NI = floor(Lx/h)
NJ = floor(Ly/h)
dx = Lx/NI
dy = Ly/NJ
h = (dx+dy)/2
NI = NI + 1       %left layer
NJ = NJ + 2       %top and bottom layer
umax = 1

%Material properties
nu = 1/100;
rho= 1;




%% Mesh definition
fprintf('generating Mesh ...\n'); %fflush(stdout);

%spacingA;
%regularMesh;

[Xp,Yp,Hxp,Hyp,NI,NJ,aux,aux2] = createMeshB ();


[ Xu,Yu,Xv,Yv,dimXu,dimYu,dimXv,dimYv ] = getXuXv (Xp,Yp,Hxp,Hyp);

plotMesh;


%% Memory allocation

%Arrays
dimP = NI*NJ
dimU = (NI+1)*NJ
dimV = NI*(NJ+1)

Ap = spalloc(dimP,dimP,5*dimP);
Dxu = spalloc(dimU,dimU,5*dimU);
Dyv = spalloc(dimV,dimV,5*dimV);

AxuE = spalloc(dimU,dimU,2*dimU);
AxuW = spalloc(dimU,dimU,2*dimU);
AxuN = spalloc(dimU,dimU,2*dimU);
AxuS = spalloc(dimU,dimU,2*dimU);
AxvN = spalloc(dimU,dimV,2*dimU);
AxvS = spalloc(dimU,dimV,2*dimU);
AxuEs = spalloc(dimU,dimU,2*dimU);
AxuWs = spalloc(dimU,dimU,2*dimU);
AxuNs = spalloc(dimU,dimU,2*dimU);
AxuSs = spalloc(dimU,dimU,2*dimU);

AyvE = spalloc(dimV,dimV,2*dimV);
AyvW = spalloc(dimV,dimV,2*dimV);
AyvN = spalloc(dimV,dimV,2*dimV);
AyvS = spalloc(dimV,dimV,2*dimV);
AyuE = spalloc(dimV,dimU,2*dimV);
AyuW = spalloc(dimV,dimU,2*dimV);
AyvEs = spalloc(dimV,dimV,2*dimV);
AyvWs = spalloc(dimV,dimV,2*dimV);
AyvNs = spalloc(dimV,dimV,2*dimV);
AyvSs = spalloc(dimV,dimV,2*dimV);

gradU = spalloc(dimP,dimU,2*dimP);
gradV = spalloc(dimP,dimV,2*dimP);
gradPU = spalloc(dimU,dimP,2*dimU);
gradPV = spalloc(dimV,dimP,2*dimV);

avgU = spalloc(dimP,dimU,2*dimP);
avgV = spalloc(dimP,dimV,2*dimP);




%% Temporal integration

hmin1 = min (Hyp);
hmin2 = min (Hxp);
hmin = min(hmin1,hmin2);
dtd = hmin*hmin/(nu*4)
dtc = 2*nu/umax
%dt = 0.001 
dt = 0.5*min(dtd,dtc)
T0 = 0;
T1 = 1000;
NT = (T1-T0)/dt;
dtol = 1.0E-6;

%Postprocessing
NPRIN = 100;
NSAVE = 0;

%% Mask definition
fprintf('generating Mask ...\n'); %fflush(stdout);

maskFernando;
%getMask;


%% Boundary Conditions
fprintf('generating boundary conditions ...\n'); %fflush(stdout);


pbc = zeros(dimP,1);
ubc = zeros(dimU,1);
vbc = zeros(dimV,1);
N=0;
for j=1:NJ
    if (masku(ij2nu(2,j)) == INFLOW)
        N = N + 1;
    end
end

for j=1:NJ
    if (masku(ij2nu(2,j)) == INFLOW)
        y = Yu(j);
        ubc(ij2nu(2,j)) = umax*(-4*y^2+12*y-8);
    end
end

% UI=parabola(N);
% N=1;
% for j=1:NJ
%     if (masku(ij2nu(2,j)) == INFLOW)
%         ubc(ij2nu(2,j)) = umax*UI(N);
%         N = N + 1;
%     end
% end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% PRESSURE LAPLACIAN
fprintf('pressure laplacian ...\n'); %fflush(stdout);


for i=1:NI
    
    for j=1:NJ
        
        % getSpacing
        [ dxr,dxl,dyt,dyb,hx,hy ] = getSpacingP (i,j,Xp,Yp,Hxp,Hyp);
        
        if (mask(ij2n(i,j)) == INTERIOR)
            %-- east face
            %fprintf('%d, %d, %d, %d \n' ,dxr,dxl,dyb,dyt)
            if (masku(ij2nu(i+1,j)) == INTERIOR)
                Ap(ij2n(i,j),ij2n(i,j)) = Ap(ij2n(i,j),ij2n(i,j)) + -1/(dxr*dxl);
                Ap(ij2n(i,j),ij2n(i+1,j)) =  dxl/(dxl*dxr*(dxl+dxr));
            end
            %-- west face
            if (masku(ij2nu(i,j)) == INTERIOR)
                Ap(ij2n(i,j),ij2n(i,j)) = Ap(ij2n(i,j),ij2n(i,j)) + -1/(dxr*dxl);
                Ap(ij2n(i,j),ij2n(i-1,j)) =  dxr/(dxl*dxr*(dxl+dxr));
            end
            %-- north face
            if (maskv(ij2nv(i,j+1)) == INTERIOR)
                Ap(ij2n(i,j),ij2n(i,j)) = Ap(ij2n(i,j),ij2n(i,j)) + -1/(dyb*dyt);
                Ap(ij2n(i,j),ij2n(i,j+1)) =  dyb/(dyb*dyt*(dyb+dyt));
            end
            %-- south face
            if (maskv(ij2nv(i,j)) == INTERIOR)
                Ap(ij2n(i,j),ij2n(i,j)) = Ap(ij2n(i,j),ij2n(i,j)) + -1/(dyb*dyt);
                Ap(ij2n(i,j),ij2n(i,j-1)) =  dyt/(dyb*dyt*(dyb+dyt));
            end
        else
            %not interior: left unchanged
            Ap(ij2n(i,j),ij2n(i,j)) = 1;
        end
    end
end





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% FACE VELOCITIES (U-Cells)
fprintf('face velocities (u-cells) ...\n'); %fflush(stdout);


for i=1:NI+1
    for j=1:NJ
        
        
        % getSpacing
        [ hxu,hyu,dvu ] = getSpacingU( i,j,Xu,Yu );
        
        % hxu, hyu : U cell size
        % dvu      : U cell volume
        

        
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%% NORMAL VELOCITIES
        %%%%   ux{E,W,N,S} = Axu{E,W,N,S}*u
        if (masku(ij2nu(i,j)) == INTERIOR)
            %-- east face
            if (masku(ij2nu(i+1,j)) == INTERIOR)
                AxuE(ij2nu(i,j),ij2nu(i,j)) = 1/2;
                AxuE(ij2nu(i,j),ij2nu(i+1,j)) = 1/2;
                AxuEs(ij2nu(i,j),ij2nu(i,j)) = 1/2 * hyu / dvu;
                AxuEs(ij2nu(i,j),ij2nu(i+1,j)) = 1/2 * hyu / dvu;
                
            elseif (masku(ij2nu(i+1,j)) == WALL)
                AxuE(ij2nu(i,j),ij2nu(i,j)) = 1/2;
                AxuEs(ij2nu(i,j),ij2nu(i,j)) = 1/2 * hyu / dvu;
            elseif (masku(ij2nu(i+1,j)) == OUTFLOW)
                % zero gradient
                AxuE(ij2nu(i,j),ij2nu(i,j)) = 1;
                AxuEs(ij2nu(i,j),ij2nu(i,j)) = 1 * hyu / dvu;
            end
            %-- west face
            if (masku(ij2nu(i-1,j)) == INTERIOR)
                AxuW(ij2nu(i,j),ij2nu(i,j)) = 1/2;
                AxuW(ij2nu(i,j),ij2nu(i-1,j)) = 1/2;
                AxuWs(ij2nu(i,j),ij2nu(i,j)) = 1/2 * hyu/ dvu;
                AxuWs(ij2nu(i,j),ij2nu(i-1,j)) = 1/2 * hyu/ dvu;
                
            elseif (masku(ij2nu(i-1,j)) == WALL)
                AxuW(ij2nu(i,j),ij2nu(i,j)) = 1/2;
                AxuWs(ij2nu(i,j),ij2nu(i,j)) = 1/2 * hyu/ dvu;
                
            elseif (masku(ij2nu(i-1,j)) == INFLOW)
                AxuW(ij2nu(i,j),ij2nu(i,j)) = 1/2;
                AxuW(ij2nu(i,j),ij2nu(i-1,j)) = 1/2;
                AxuWs(ij2nu(i,j),ij2nu(i,j)) = 1/2 * hyu/ dvu;
                AxuWs(ij2nu(i,j),ij2nu(i-1,j)) = 1/2 * hyu/ dvu;
            end
            %--north face
            if (masku(ij2nu(i,j+1)) == INTERIOR)
                AxuN(ij2nu(i,j),ij2nu(i,j)) = 1/2;
                AxuN(ij2nu(i,j),ij2nu(i,j+1)) = 1/2;
                AxuNs(ij2nu(i,j),ij2nu(i,j)) = 1/2 * hxu/ dvu;
                AxuNs(ij2nu(i,j),ij2nu(i,j+1)) = 1/2 * hxu/ dvu;
            elseif (masku(ij2nu(i,j+1)) == WALL)
                %do nothing
            end
            %--south face
            if (masku(ij2nu(i,j-1)) == INTERIOR)
                AxuS(ij2nu(i,j),ij2nu(i,j)) = 1/2;
                AxuS(ij2nu(i,j),ij2nu(i,j-1)) = 1/2;
                AxuSs(ij2nu(i,j),ij2nu(i,j)) = 1/2 * hxu/ dvu;
                AxuSs(ij2nu(i,j),ij2nu(i,j-1)) = 1/2 * hxu/ dvu;
            elseif (masku(ij2nu(i,j-1)) == WALL)
                %do nothing
            end
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %%%% TANGENCIAL VELOCITIES
            %%%%   vx{N,S}     = Axv{N,S}*v
            %--north face
            if (masku(ij2nu(i,j+1)) == INTERIOR)
                AxvN(ij2nu(i,j),ij2nv(i,j+1)) = 1/2;
                AxvN(ij2nu(i,j),ij2nv(i-1,j+1)) = 1/2;
            elseif (masku(ij2nu(i,j+1)) == WALL)
                %do nothing
            end
            %--south face
            if (masku(ij2nu(i,j-1)) == INTERIOR)
                AxvS(ij2nu(i,j),ij2nv(i-1,j)) = 1/2;
                AxvS(ij2nu(i,j),ij2nv(i,j)) = 1/2;
            elseif (masku(ij2nu(i,j-1)) == WALL)
                %do nothing
            end
        end
        
    end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% FACE VELOCITIES (V-Cells)
fprintf('face velocities (v-cells) ...\n'); %fflush(stdout);


for i=1:NI
    for j=1:NJ+1
        
        
        % getSpacing
        [ hxv,hyv,dvv ] = getSpacingV( i,j,Xv,Yv );
        
        
        
        %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
        %%%% NORMAL VELOCITIES
        %%%%   vy{E,W,N,S} = Ayv{E,W,N,S}*v
        if (maskv(ij2nv(i,j)) == INTERIOR)
            %--east face
            if (maskv(ij2nv(i+1,j)) == INTERIOR)
                AyvE(ij2nv(i,j),ij2nv(i,j)) = 1/2;
                AyvE(ij2nv(i,j),ij2nv(i+1,j)) = 1/2;
                AyvEs(ij2nv(i,j),ij2nv(i,j)) = 1/2 * hyv / dvv;
                AyvEs(ij2nv(i,j),ij2nv(i+1,j)) = 1/2 * hyv / dvv;
            elseif (maskv(ij2nv(i+1,j)) == WALL)
                %do nothing
            elseif (maskv(ij2nv(i+1,j)) == OUTFLOW)
                AyvE(ij2nv(i,j),ij2nv(i,j)) = 1/2;
                AyvEs(ij2nv(i,j),ij2nv(i,j)) = 1/2 * hyv / dvv;
            end
            %--west face
            if (maskv(ij2nv(i-1,j)) == INTERIOR)
                AyvW(ij2nv(i,j),ij2nv(i,j)) = 1/2;
                AyvW(ij2nv(i,j),ij2nv(i-1,j)) = 1/2;
                AyvWs(ij2nv(i,j),ij2nv(i,j)) = 1/2 * hyv / dvv;
                AyvWs(ij2nv(i,j),ij2nv(i-1,j)) = 1/2 * hyv / dvv;
            elseif (maskv(ij2nv(i-1,j)) == WALL)
                %do nothing
            elseif (maskv(ij2nv(i-1,j)) == INFLOW)
                %do nothing
            end
            %--north face
            if (maskv(ij2nv(i,j+1)) == INTERIOR)
                AyvN(ij2nv(i,j),ij2nv(i,j)) = 1/2;
                AyvN(ij2nv(i,j),ij2nv(i,j+1)) = 1/2;
                AyvNs(ij2nv(i,j),ij2nv(i,j)) = 1/2 * hxv / dvv;
                AyvNs(ij2nv(i,j),ij2nv(i,j+1)) = 1/2 * hxv / dvv;
            elseif (maskv(ij2nv(i,j+1)) == WALL)
                AyvN(ij2nv(i,j),ij2nv(i,j)) = 1/2;
                AyvNs(ij2nv(i,j),ij2nv(i,j)) = 1/2 * hxv / dvv;
            end
            %--south face
            if (maskv(ij2nv(i,j-1)) == INTERIOR)
                AyvS(ij2nv(i,j),ij2nv(i,j)) = 1/2;
                AyvS(ij2nv(i,j),ij2nv(i,j-1)) = 1/2;
                AyvSs(ij2nv(i,j),ij2nv(i,j)) = 1/2 * hxv / dvv;
                AyvSs(ij2nv(i,j),ij2nv(i,j-1)) = 1/2 * hxv / dvv;
            elseif (maskv(ij2nv(i,j-1)) == WALL)
            end
            %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
            %%%% TANGENCIAL VELOCITIES
            %%%%   uy{E,W}     = Ayu{E,W}*u
            %--east face
            if (maskv(ij2nv(i+1,j)) == INTERIOR)
                AyuE(ij2nv(i,j),ij2nu(i+1,j)) = 1/2;
                AyuE(ij2nv(i,j),ij2nu(i+1,j-1)) = 1/2;
            elseif (maskv(ij2nv(i+1,j)) == WALL)
                %do nothing
            elseif (maskv(ij2nv(i+1,j)) == OUTFLOW)
                AyuE(ij2nv(i,j),ij2nu(i+1,j)) = 1/2;
                AyuE(ij2nv(i,j),ij2nu(i+1,j-1)) = 1/2;
            end
            %--west face
            if (maskv(ij2nv(i-1,j)) == INTERIOR)
                AyuW(ij2nv(i,j),ij2nu(i,j)) = 1/2;
                AyuW(ij2nv(i,j),ij2nu(i,j-1)) = 1/2;
            elseif (maskv(ij2nv(i-1,j)) == WALL)
                %do nothing
            elseif (maskv(ij2nv(i-1,j)) == INFLOW)
                AyuW(ij2nv(i,j),ij2nu(i,j)) = 1/2;
                AyuW(ij2nv(i,j),ij2nu(i,j-1)) = 1/2;
            end
        end
    end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% VISCOUS LAPLACIAN (U-Cells)
fprintf('viscous laplacian (u-cells) ...\n'); %fflush(stdout);



for i=1:NI+1
    
    for j=1:NJ
        
        [ hxu,hyu,dvu,dxur,dxul,dyut,dyub ] = getSpacingU( i,j,Xu,Yu );
        
        
        
        if (masku(ij2nu(i,j)) == INTERIOR)
            %-- east face
            if (masku(ij2nu(i+1,j)) == INTERIOR)
                Dxu(ij2nu(i,j),ij2nu(i,j)) =  Dxu(ij2nu(i,j),ij2nu(i,j)) + -hyu/(dvu*dxur);
                Dxu(ij2nu(i,j),ij2nu(i+1,j)) =  hyu/(dvu*dxur);
            elseif (masku(ij2nu(i+1,j)) == WALL)
                Dxu(ij2nu(i,j),ij2nu(i,j)) =  Dxu(ij2nu(i,j),ij2nu(i,j)) + -hyu/(dvu*dxur);
            elseif (masku(ij2nu(i+1,j)) == OUTFLOW)
                %zero gradient
                %do nothing
            end
            %-- west face
            if (masku(ij2nu(i-1,j)) == INTERIOR)
                Dxu(ij2nu(i,j),ij2nu(i,j)) =  Dxu(ij2nu(i,j),ij2nu(i,j)) + -hyu/(dvu*dxul);
                Dxu(ij2nu(i,j),ij2nu(i-1,j)) =  hyu/(dvu*dxul);
            elseif (masku(ij2nu(i-1,j)) == WALL)
                Dxu(ij2nu(i,j),ij2nu(i,j)) =  Dxu(ij2nu(i,j),ij2nu(i,j)) + -hyu/(dvu*dxul);
            elseif (masku(ij2nu(i-1,j)) == INFLOW)
                Dxu(ij2nu(i,j),ij2nu(i,j)) =  Dxu(ij2nu(i,j),ij2nu(i,j)) + -hyu/(dvu*dxul);
                Dxu(ij2nu(i,j),ij2nu(i-1,j)) =  hyu/(dvu*dxul);
            end
            %-- north face
            if (masku(ij2nu(i,j+1)) == INTERIOR)
                Dxu(ij2nu(i,j),ij2nu(i,j)) =  Dxu(ij2nu(i,j),ij2nu(i,j)) + -hxu/(dvu*dyut);
                Dxu(ij2nu(i,j),ij2nu(i,j+1)) =  hxu/(dvu*dyut);
            elseif (masku(ij2nu(i,j+1)) == WALL)
                %ghost point
                Dxu(ij2nu(i,j),ij2nu(i,j)) =  Dxu(ij2nu(i,j),ij2nu(i,j)) + -2*hxu/(dvu*dyut);
            end
            %-- south face
            if (masku(ij2nu(i,j-1)) == INTERIOR)
                Dxu(ij2nu(i,j),ij2nu(i,j)) =  Dxu(ij2nu(i,j),ij2nu(i,j)) + -hxu/(dvu*dyub);
                Dxu(ij2nu(i,j),ij2nu(i,j-1)) =  hxu/(dvu*dyub);
            elseif (masku(ij2nu(i,j-1)) == WALL)
                %ghost point
                Dxu(ij2nu(i,j),ij2nu(i,j)) =  Dxu(ij2nu(i,j),ij2nu(i,j)) + -2*hxu/(dvu*dyub);
            end
        else
            %not interior: left unchanged
            Dxu(ij2nu(i,j),ij2nu(i,j)) = 1;
        end
    end
end




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% VISCOUS LAPLACIAN (V-Cells)
fprintf('viscous laplacian (v-cells) ...\n'); %fflush(stdout);



for i=1:NI
    
    for j=1:NJ+1
        
        
        
        [ hxv,hyv,dvv,dxvr,dxvl,dyvt,dyvb ] = getSpacingV( i,j,Xv,Yv );
        
        
        
        if (maskv(ij2nv(i,j)) == INTERIOR)
            %-- east face
            if (maskv(ij2nv(i+1,j)) == INTERIOR)
                Dyv(ij2nv(i,j),ij2nv(i,j))  =  Dyv(ij2nv(i,j),ij2nv(i,j)) + -hyv/(dvv*dxvr);
                Dyv(ij2nv(i,j),ij2nv(i+1,j)) =  hyv/(dvv*dxvr);
            elseif (maskv(ij2nv(i+1,j)) == WALL)
                %ghost point
                Dyv(ij2nv(i,j),ij2nv(i,j))  =  Dyv(ij2nv(i,j),ij2nv(i,j)) + -2*hyv/(dvv*dxvr);
            elseif (maskv(ij2nv(i+1,j)) == OUTFLOW)
                Dyv(ij2nv(i,j),ij2nv(i,j))  =  Dyv(ij2nv(i,j),ij2nv(i,j)) + -hyv/(dvv*dxvr);
            end
            %-- west face
            if (maskv(ij2nv(i-1,j)) == INTERIOR)
                Dyv(ij2nv(i,j),ij2nv(i,j))  =  Dyv(ij2nv(i,j),ij2nv(i,j)) + -hyv/(dvv*dxvl);
                Dyv(ij2nv(i,j),ij2nv(i-1,j)) =  hyv/(dvv*dxvl);
            elseif (maskv(ij2nv(i-1,j)) == WALL)
                %ghost point
                Dyv(ij2nv(i,j),ij2nv(i,j))  =  Dyv(ij2nv(i,j),ij2nv(i,j)) + -2*hyv/(dvv*dxvl);
            elseif (maskv(ij2nv(i-1,j)) == INFLOW)
                %ghost point
                Dyv(ij2nv(i,j),ij2nv(i,j))  =  Dyv(ij2nv(i,j),ij2nv(i,j)) + -2*hyv/(dvv*dxvl);
            end
            %-- north face
            if (maskv(ij2nv(i,j+1)) == INTERIOR)
                Dyv(ij2nv(i,j),ij2nv(i,j))  =  Dyv(ij2nv(i,j),ij2nv(i,j)) + -hxv/(dvv*dyvt);
                Dyv(ij2nv(i,j),ij2nv(i,j+1)) =  hxv/(dvv*dyvt);
            elseif (maskv(ij2nv(i+1,j)) == WALL)
                Dyv(ij2nv(i,j),ij2nv(i,j))  =  Dyv(ij2nv(i,j),ij2nv(i,j)) + -hxv/(dvv*dyvt);
            end
            %-- south face
            if (maskv(ij2nv(i,j-1)) == INTERIOR)
                Dyv(ij2nv(i,j),ij2nv(i,j))  =  Dyv(ij2nv(i,j),ij2nv(i,j)) + -hxv/(dvv*dyvb);
                Dyv(ij2nv(i,j),ij2nv(i,j-1)) =  hxv/(dvv*dyvb);
            elseif (maskv(ij2nv(i+1,j)) == WALL)
                Dyv(ij2nv(i,j),ij2nv(i,j))  =  Dyv(ij2nv(i,j),ij2nv(i,j)) + -hxv/(dvv*dyvb);
            end
        else
            %not interior: left unchanged
            Dyv(ij2nv(i,j),ij2nv(i,j)) = 1;
        end
    end
end




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% VELOCITY GRADIENTS (P-Cells)
fprintf('velocity gradients ...\n'); %fflush(stdout);


for i=1:NI
    for j=1:NJ
        
        % getSpacing
        [ dxr,dxl,dyt,dyb,hx,hy ] = getSpacingP (i,j,Xp,Yp,Hxp,Hyp);
        
        if (mask(ij2n(i,j)) == INTERIOR)
            gradU(ij2n(i,j),ij2nu(i+1,j)) = 1/hx;
            gradU(ij2n(i,j),ij2nu(i,j)) = -1/hx;
            gradV(ij2n(i,j),ij2nv(i,j+1)) = 1/hy;
            gradV(ij2n(i,j),ij2nv(i,j)) = -1/hy;
        end
    end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% PRESSURE GRADIENT (U-Cells)
fprintf('pressure gradient (u-cells) ...\n'); %fflush(stdout);

for i=1:NI+1
    for j=1:NJ
        
        [ hxu,hyu,dvu,dxur,dxul,dyut,dyub ] = getSpacingU( i,j,Xu,Yu );
        
        if (masku(ij2nu(i,j)) == INTERIOR)
            gradPU(ij2nu(i,j),ij2n(i,j)) = 1/hxu;
            gradPU(ij2nu(i,j),ij2n(i-1,j)) = -1/hxu;
        end
    end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% PRESSURE GRADIENT (V-Cells)
fprintf('pressure gradient (v-cells) ...\n'); %fflush(stdout);


for i=1:NI
    for j=1:NJ+1
        
        [ hxv,hyv,dvv,dxvr,dxvl,dyvt,dyvb ] = getSpacingV( i,j,Xv,Yv );
        
        if (maskv(ij2nv(i,j)) == INTERIOR)
            gradPV(ij2nv(i,j),ij2n(i,j)) = 1/hyv;
            gradPV(ij2nv(i,j),ij2n(i,j-1)) = -1/hyv;
        end
    end
end


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% VELOCITY AVERAGES (for visualization)
fprintf('velocity averages ...\n'); %fflush(stdout);


for i=1:NI
    for j=1:NJ
        if (mask(ij2n(i,j)) == INTERIOR)
            avgU(ij2n(i,j),ij2nu(i,j)) = 1/2;
            avgU(ij2n(i,j),ij2nu(i+1,j)) = 1/2;
            avgV(ij2n(i,j),ij2nv(i,j)) = 1/2;
            avgV(ij2n(i,j),ij2nv(i,j+1)) = 1/2;
        end
    end
end

%%% END INITIALIZATION
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% POST-PROCESSING INITIALIZATION
XP=Xp;
YP=Yp;
XU=Xu;
YU=Yu;

XV=Xv;
YV=Yv;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% BEGIN TIME INTEGRATION

p = zeros(dimP,1);
u = zeros(dimU,1);
v = zeros(dimV,1);

u = ubc;

i = 0;
t = T0;
while ((t < T1)&&(i < NT))
    
    %advection operator
    uEs = AxuEs*u;
    uWs = AxuWs*u;
    uNs = AxuNs*u;
    uSs = AxuSs*u;
    
    uE = AxuE*u; uW = AxuW*u;
    uN = AxuN*u; uS = AxuS*u;
    vN = AxvN*v; vS = AxvS*v;
    Ax = (uEs.*uE - uWs.*uW + uNs.*vN - uSs.*vS);
    
    vNs = AyvNs*v;
    vSs = AyvSs*v;
    vEs = AyvEs*v;
    vWs = AyvWs*v;
    
    vN = AyvN*v;
    vS = AyvS*v;
    vE = AyvE*v;
    vW = AyvW*v;
    uE = AyuE*u;
    uW = AyuW*u;
    Ay = (vNs.*vN - vSs.*vS + vEs.*uE - vWs.*uW);
    
    %viscous operator
    Dx = Dxu*u-ubc;
    Dy = Dyv*v-vbc;
    
    %prediction step
    ustar = u + dt*(-Ax + nu*Dx);
    vstar = v + dt*(-Ay + nu*Dy);
    
    %pressure poisson equation
    rhs = (rho/dt)*(gradU*ustar + gradV*vstar);
    pnew = Ap\rhs;
    
    %correction step
    unew = ustar - (dt/rho)*gradPU*pnew;
    vnew = vstar - (dt/rho)*gradPV*pnew;
    
    %deltas
    dp=norm(pnew-p);
    du=norm(unew-u);
    dv=norm(vnew-v);
    
    %finish timestep
    u = unew;
    v = vnew;
    p = pnew;
    t = t + dt;
    i = i + 1;
    
    %output ?
    if (mod(i,NPRIN) == 0)
        fprintf('t= %12.5E  |u|= %12.5E |v|= %12.5E |p|= %12.5E\n',...
            t,norm(u),norm(v),norm(p));
        
        %fflush(stdout);
    end
    
    %save image ?
    if (mod(i,NSAVE) == 0)
        fprintf('saving solution at T= %g\n', t);
        %fflush(stdout);
        
        P=reshape(p,NJ,NI);
        U=reshape(u,NJ,NI+1);
        V=reshape(v,NJ+1,NI);
        
        figure(1,'visible','off');
        clf();
        contourf(XP(2:NI),YP(2:NJ-1),P(2:NJ-1,2:NI));
        axis([XU(1) XU(NI+1) YV(1) YV(NJ+1)])
        colorbar();
        title(sprintf('Pressure at T=%g',t));
        print(sprintf('press%08d.png',i),'-dpng');
        close();
        
        figure(2,'visible','off');
        clf();
        contourf(XU(2:NI),YU(2:NJ-1),U(2:NJ-1,2:NI));
        hold
        czero=contour(XU(2:NI),YU(2:NJ-1),U(2:NJ-1,2:NI),[0,0]);
        axis([XU(1) XU(NI+1) YV(1) YV(NJ+1)])
        colorbar();
        xzero=czero(end-1)-dx;
        title(sprintf('U-Velocity at T=%g, Vortex Attachment Point at X=%g'),...
            t,xzero);
        print(sprintf('velou%08d.png',i),'-dpng');
        close();
        
        figure(3,'visible','off');
        clf();
        contourf(XV(2:NI),YV(2:NJ),V(2:NJ,2:NI));
        axis([XU(1) XU(NI+1) YV(1) YV(NJ+1)])
        colorbar();
        title(sprintf('V-Velocity at T=%g',t));
        print(sprintf('velov%08d.png',i),'-dpng');
        close();
    end
    
    %convergence ?
    if ((dp < dtol)&&(du < dtol)&&(dv < dtol))
        break;
    end
end



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% POST-PROCESSING
P=reshape(p,NJ,NI);
U=reshape(u,NJ,NI+1);
V=reshape(v,NJ+1,NI);

% for i=1:NI+1
%     U(NJ,i) = 0.5*(U (NJ,1)+U(NJ-1,i));
%     YU (NJ) = 0.5*(Yu (NJ)+Yu(NJ-1));
% end
% 
% for i=1:NI+1
%     U(1,i) = 0.5*(U (2,1)+U(2,i));
%     YU (1) = 0.5*(Yu (2)+Yu(2));
% end



figure();
contour(XP(2:NI),YP(2:NJ-1),P(2:NJ-1,2:NI));
axis([XU(1) XU(NI+1) YV(1) YV(NJ+1)])
colorbar();
title(sprintf('Pressure at T=%g',t));

figure();

contourf(XU(1:NI),YU(1:NJ),U(1:NJ,1:NI));
hold
pbaspect ([12 2 1])
czero=contour(XU(2:NI),YU(2:NJ),U(2:NJ,2:NI),[0,0]);
axis([XU(1) XU(NI+1) YV(1) YV(NJ+1)])
colorbar();
xzero=czero(end-1)-dx;
title(sprintf('U-Velocity at T=%g, Vortex Attachment Point at X=%g',...
    t,xzero));

figure();
contourf(XV(2:NI),YV(2:NJ),V(2:NJ,2:NI));
axis([XU(1) XU(NI+1) YV(1) YV(NJ+1)])
colorbar();
title(sprintf('V-Velocity at T=%g',t));

figure();
up = avgU*u;
vp = avgV*v;
UP=reshape(up,NJ,NI);
VP=reshape(vp,NJ,NI);

h=quiver(XP(2:NI),YP(2:NJ-1),UP(2:NJ-1,2:NI),VP(2:NJ-1,2:NI));
pbaspect ([12 2 1])
axis([XU(1) XU(NI+1) YV(1) YV(NJ+1)])
set (h, 'maxheadsize', 0.33);

%%% END POST-PROCESSING
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
